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Abstract 

Orientation and microstructure are investigated in sheared Brownian suspensions of 
hard dicolloidal particles, with the dicolloids modeled as two fused spheres of varying 
radii and center to center separations. Two different particle shapes named homonuclear 
(aspect ratio 1.1) and fused-dumbbells (aspect ratio 1.5) were considered. Hydrodynamic 
interactions between the particles were computed with a modified lubrication model called 
Fast Lubrication Dynamics. Studies were conducted for a wide range of volume fractions 
between 0.3 < (j) < 0.5 and Peclet numbers between < Pe < 1000. The microstructure 
was found to be disordered at all volume fractions, though signatures of weak string like 
ordering were evident particularly in — 0.5 homonuclear suspensions at intermediate to 
high shear rates (Pe in the range 10 — 100). Complex orientation behavior was observed 
as a function of shape, shear rates, and volume fractions. At very low shear rates, random 
orientation distribution was observed in all cases. At the highest shear rates, orientation 
distribution in suspensions of homonuclear particles exhibited a shift towards an alignment 
with the vorticity axis at all volume fractions, while in suspensions of fused-dumbbells it 
exhibited a shift away from the vorticity axis at low volume fractions and a negligible 
shift at higher volume fractions. The orientation behavior is further characterized by 
examining the orientation distribution in the velocity-gradient plane - in this case an 
increased particle alignment with the velocity axis is generally observed with increasing 
volume fractions, but not universally with increasing shear rates. Mechanistic origins for 
the complex orientation behavior as a function of shear rate, volume fraction, and particle 
shape is described. 

1 Introduction 

In recent times anisotropic particles have attracted considerable interest in a wide range of 
fields. This surge in interest is arguably fueled by the recent synthesis of a variety of novel 
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Figure 1: Examples of dicolloidal particles investigated in the current work. The two particle 
shapes are named homonuclear and fused-dumbbell particles as noted in the figure. 

particle shapes which themselves are enabled by the recent advances in particle synthesis 
techniques. Sizes of these particles are usually in the nanometer range with O(IOO) nm being 
a typical size. Many applications have been proposed for these anisotropic particles; some 
prominent examples include their assembly into advanced nanostructured materials'^ and for 
drug delivery ^'^ among others. An interesting class of recently synthesized anisotropic particles 
are the dicolloidal particles whose geometry is closely approximated by the union of two 
intersecting spheres of varying radii and center to center separations. Figure (1) illustrates 
several dicolloidal particle shapes examined in the present effort. A significant feature of 
dicolloidal particles is its low aspect ratio (AR) with an upper bound of 2.0. Owing to this low 
aspect ratio, mildly aspherical dicolloids often exhibit microstructural and rheological behavior 
similar to that in suspensions of spheres. However, in addition to the positional microstructure, 
suspensions of dicolloids possess distinct orientational microstructure. The goals of present 
work are to characterize both the orientational and the positional microstructure in sheared 
suspensions of Brownian dicolloidal particles. This continues our previous work on dicolloidal 
particles where we investigated the orientational and positional microstructure in sheared non- 
Brownian suspensions^. We review some of the important observations in this study below. 
In addition, we review below other pertinent literature on the orientation behavior in sheared 
suspensions of anisotropic particles. 

In the absence of non-hydrodynamic effects (such as Brownian motion or electrostatic 
forces), Jeffery^*^ showed that a single spheroid (or more generally an ellipsoid) in a linear 
shear flow would rotate indefinitely in a single parameter family of closed orbits, commonly 
known as the Jeffery orbits. The parameter characterizing an orbit is called the orbit constant 
C whose value varies from C = for particle director aligned with the vorticity axis to C = oo 
for particle alignment in the velocity-gradient plane. For an axisymmetric body, such as a 
spheroid or a dicolloid, its orientation may be specified by spherical coordinates (if, 9) where 
measures the angle between the particle director and z axis, while (p is angle between the 
director and the y axis measured clockwise in the x — y plane (see Fig. 2). The evolution of 
the orientational angles ((/?, 9) with time t is given in terms of the orbit constant C by 
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Figure 2: Definition of the spherical polar coordinate 9 and the azimuthal angle 
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where 7 is the shear rate, is the aspect ratio of the spheroid, while k denotes the initial 
phase in the orbit. A particularly simple limit is that of a spherical particle (r^ = 1), in which 
case Jeffery's orbits corresponds to lines of constant 6 (more specifically, 9 = tan~^ C). Later, 
Bretherton^^ showed that in a linear shear flow most axisymmetric bodies, which includes 
dicolloids, will also rotate in JefFery orbits with in the equations above replaced by an 
effective aspect ratio r^. of the particle. It can be shown that the rate of rotation of a particle 
in any given orbit is non-uniform, such that its director spends a greater fraction of the time 
aligned in the velocity-vorticity plane in comparison to the gradient-vorticity plane. As such, 
under shear, a single non-Brownian anisotropic particle by itself will exhibit an enhanced 
time-averaged flow alignment (parallel to the velocity axis) in conjunction with a reduced 
gradient alignment (henceforth, by particle alignment we imply particle director alignment). 
On the other hand, the time-averaged vorticity alignment is mostly determined by the initial 
orientation of the particle; the initial orientation sets the orbit constant C, and the vorticity 
alignment does not change appreciably with time in any given orbit for particles with small 
degree of anisotropy. 

The effect of Brownian motion on the orientation behavior of axisymmetric particles have 
been characterized in several theoretical studies in the limit of infinite dilution, i.e., non- 
interacting particle suspensions. Leal and Hinch^^ investigated the orientation distribution 
in suspensions of spheroids in the limit of very weak Brownian motion. The key argument 
in this study was that the rotation of a particle in any given Jeffery's orbit is unaltered 
when the Brownian motion is very weak. However, after sufficiently long times, the first 
effect of the weak Brownian motion will be to yield steady state distributions across orbits 
independent of the initial state of the suspension. Since the rate of rotation of a particle 
in any given orbit is unchanged, this result implies an enhanced particle alignment with the 
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velocity axis and a reduced particle alignment with the gradient axis, say, in comparison to 
an isotropic distribution. In an earlier study^, we showed that in the limit of weak Brownian 
motion, the particle alignment with the vorticity axis remains approximately isotropic for 
dicolloids. Hinch and Leal^'^ later extended their previous analysis to include higher order 
approximations, though still in the weak Brownian limit. This study, in comparison to Jeffery's 
solution, showed an enhanced probability at orientations just prior to complete flow alignment, 
in conjunction with a reduced probability at orientations just after complete flow alignment. 
This behavior results from the fact that the Brownian motion causes diffusion away from 
the highly populated flow aligned states, which, to first approximation, occurs with equal 
probability towards the pre flow aligned state and post flow aligned state. The advection, on 
the other hand, causes particle motion in a single direction of the undisturbed Jeffery's orbit. 
The balance between these two effects leads to the aforementioned probability distribution. 

In the limit of strong Brownian motion, the orientation distribution in non-interacting 
suspensions of axisymmetric particles is provided by, e.g., Brenner At equilibrium, as 
could be expected, the orientation distribution was found to be isotropic. The first effect of 
the weak flow is to enhance the orientation probability around the extensional axis of the 
shear flow and a reduced probability along its compressional axis. The effect of shear flow 
is understood by noting that the overall flow in this case can be expressed as a sum of an 
extensional component and a rotational component. The rotational component has no effect 
at the leading order as the rate of rotation of the particle due to this component is uniform and 
the initial distribution, i.e., at equilibrium, isotropic. The effect of the extensional component 
of the flow is to rotate the particle away from the compressional axis towards alignment with 
the extensional axis, which explains the above observation. At higher order approximations, 
both the extensional as well as the rotational component of the flow will contribute. The effect 
of the rotational component of the flow will be to rotate the orientation distribution clockwise, 
thereby bringing the probability maxima closer to the velocity axis and the probability minima 
closer to the gradient axis. At very high shear rates, it can be expected that the weak Brownian 
limit discussed above will be reached, i.e., the region of maximum orientation probability 
will coincide with the velocity axis, while the region of minimum orientation probability will 
coincide with the gradient axis. 

There have been relatively few studies on the orientation behavior in Brownian suspensions 
of anisotropic particles at finite concentrations. Some are discussed next. Meng and Higdon^^ 
investigated the orientation behavior in concentrated Brownian suspensions of plate-like par- 
ticles {AR between 3 and 7) via an extension of the Stokesian Dynamics technique. Their 
studies indicate an increasing particle alignment in the velocity- vorticity plane with increasing 
shear rates, volume fractions and aspect ratios. Egres et al. experimentally measured flow 
alignment in suspensions of ellipsoids of comparatively lower aspect ratios {AR between 2 and 
7). This study showed that the flow alignment (parallel to velocity axis) was preferred at all 
shear rates with the maximum in flow alignment coinciding with the onset of shear thickening. 



4 



At even higher shear rates, a decrease in flow ahgnment was observed. This study also indi- 
cated an increased flow alignment with higher particle loadings as well as aspect ratios. The 
authors attributed the decrease in flow alignment at high shear rates to the formation of hydro- 
clusters, which will have a smaller degree of anisotropy than the constituent particle itself; this 
reduction in the effective anisotropy explains the decrease in particle flow alignment. Besides 
the studies noted above, there have also been many other studies on the orientation behavior 
in suspensions of high aspect ratio non-Brownian fibers, which is of limited interest here. For 
a detailed survey of such suspensions the reader is referred to our earlier publication^. 

We recently investigated the orientation behavior in sheared non-Brownian suspensions 
of dicolloidal particles using the Stokesian dynamics technique^. This study explored the 
effect of particle shape as well as volume fraction. Two of the particle shapes investigated in 
this study are shown in Fig. (1). As noted in the figure, the two particles shapes, both of 
which are symmetric dicolloids, are named homonuclear and fused-dumbbell particles. The 
geometric aspect ratio of these two particles are AR = 1.1 and AR = 1.5, respectively. In this 
this study it was shown that, at low volume fractions, both the particle shapes exhibited an 
orientation drift away from the vorticity axis towards the velocity-gradient plane. This drift 
was considerably higher for the fused-dumbbell particle than the homonuclear particle. The 
origin of this behavior was attributed to non-uniform fluctuations in the orientation space 
- particles aligned with the vorticity axis were found to undergo stronger fluctuations in 
its orientation as it underwent collisions with other particles. As discussed in Shaqfeh and 
Koch^^, non-uniform fluctuations in the orientation space will result in a drift away from 
orientations with stronger fluctuations in agreement with the observations. At higher volume 
fractions, an orientation drift towards the vorticity axis away from the velocity-gradient plane 
was observed. The drift towards the vorticity axis was higher for the homonuclear particles 
than for the fused-dumbbell particles. In suspensions of force free and torque free particles, 
as is approximately the case in non-Brownian suspensions with a very short range repulsive 
force, the particle angular velocity can be directly related to the hydrodynamic stresslets, 
which refers to the symmetric part of the first moment of the fluid force on the particle. 
Indeed, it was shown that the particle drift towards the vorticity axis was correlated with 
the hydrodynamic stresslet, in particular with the second normal difference component of the 
stresslet^. In non-Brownian suspensions, the second normal stress differences are negative 
which was shown to drive the orientation distribution towards an enhanced vorticity alignment. 
It is interesting to note that, in a non-Newtonian second order fluid, Leal^^ showed that a 
high aspect ratio fiber subjected to shear will also drift towards vorticity alignment if second 
normal stress coefficient of the fluid is negative. 

In the present effort, we will conduct numerical simulations for the dynamics of sheared 
suspensions of Brownian dicolloids. The hydrodynamic interactions as well as the Brown- 
ian motion of the particle are computed using the Fast Lubrication Dynamics (FLD) algo- 
rithm The FLD algorithm is an approximation of the well known Stokesian Dynamics 
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technique ' . In the FLD algorithm, the overah resistance tensor is expressed as the sum of a 
pairwise lubrication term and a novel isotropic resistance term chosen so as to match the mean 
particle mobility from the more detailed Stokesian Dynamics. Using this algorithm we inves- 
tigate the orientation and microstructure in Brownian hard dicolloidal particle suspensions 
as a function of shear rate and volume fraction. Two different particle shapes are considered 
here: the homonuclear and the fused-dumbbell particle. The microstructure is found to lack 
any long range order in the range of volume fraction investigated (0.3 < (p < 0.5), though 
hints of weak string like ordering becomes evident at the highest volume fraction over a range 
of shear rates, especially in suspensions of homonuclear particles. The particle orientation is 
found to be random at low shear rates in all cases. At the highest shear rates, orientation 
distribution in suspensions of homonuclear particles exhibits a shift towards an alignment with 
the vorticity axis at all volume fractions, while in suspensions of fused-dumbbells it exhibits a 
shift away from the vorticity axis at low volume. The orientation distribution in the velocity- 
gradient plane shows an increased flow alignment with increasing volume fractions for both 
the particle shapes. The flow alignment is also found to increase with increasing shear rates 
at lower volume fractions, though, at higher volume fractions, this is not necessarily the case. 
The mechanism responsible for these orientation distributions with shear rate, volume frac- 
tion, and particle shape is described in light of previous studies on non-Brownian and strongly 
Brownian suspensions. 

The organization of this article is as follows. In Sec. (2) we present the Fast Lubrication 
Dynamics technique for computing the hydrodynamic and Brownian forces on the particles. 
The solution procedure for computing the unknown velocities and angular velocities of the 
particles is also discussed here. Next, in Sec. (3.1), the results for the microstructure as a 
function of volume fraction, shear rate and particle shape is presented. Following this, the 
orientation behavior is characterized in Sec. (3.2). Concluding remarks are lastly presented 
in Sec. (4). 

2 Formulation 

The motion of Brownian particles in a viscous fluid can be described by Newton's second law 
written as 

m.^=F^ + F^ + F^, (2) 

where m is the mass/moment of inertia tensor of the particles, U is the generalized veloc- 
ity/angular velocity vector, while F^, F^, and F^ are the generalized force/torque vectors 
arising respectively from the hydrodynamic stress in the fluid, stochastic Brownian motion, 
and interparticle interactions. Here the vectors U, F^, F-^ and F^ are 6A'p vectors where Np 
is the number of particles in the system; the mass/moment of inertia m is a 6Np square block 
diagonal matrix where each 6x6 block consists of the mass and moment of inertia tensor of 
each individual particle. We assume that the particle size is small, and the inertial effects are 
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negligible. With this assumption, the left hand side of (2) is zero, and the sum of the forces 
and torques on each individual particle must be zero at every instant of time: 

+ + = 0. (3) 

We discuss next the numerical calculation of the above three type of forces, i.e., F^, F^, and 
F^. 



2.1 Hydrodynamic force (F^) 

Owing to its small size, the fluid motion around the particles is governed by the Stokes 
equations. In these low Reynolds number flows, the hydrodynamic force on each particle is a 
linear function of the fluid velocity and may be determined from the solution of the governing 
equations for a specified set of particle configurations, particle velocities and of a prescribed 
undisturbed fiow field. The linear relationship may be expressed in terms of a configuration 
dependent Np body resistance tensor R as follows: 



R- • (4) 





In the above equation, a homogeneous undisturbed linear shear flow is assumed which deter- 
mines the generalized velocity /angular velocity U°° evaluated at the center of each particle 
and the uniform rate of strain tensor denoted by E°°. F^ is the usual force/torque 6iVp vector 
and denotes the stresslets on the particle, the symmetric part of the first moment of the 
force on each particle, and can be expressed by a 5Np vector. 

In a previous effort we developed a Stokesian dynamics technique to approximate the 
resistance tensor R in suspensions of non-spherical particles including dicolloids. In the Stoke- 
sian Dynamics approach, the total resistance tensor R in Eq. (4) is expressed by the following 
approximation: 

R = R^^^ + R^^, (5) 

where R^-'^^ is a many-body resistance tensor accurate for widely separated particles. When 
particles come near contact however, the small particle gaps lead to strong interactions with R 
diverging at contact. The many body resistance tensor cannot capture this singular behavior 
due to the truncated multipole expansion employed in its calculation, hence, a correction term 
is added to account for the missing terms. This correction tensor denoted by R-^^ is based 
on asymptotic lubrication theory for nearly touching particles. Unlike R*^^, R^^ is a sparse 
matrix as the lubrication interactions affect only near neighbors, and most importantly, the 
asymptotic lubrication contributions are pairwise additive as the interaction is highly localized 
around the point of contact. In the Stokesian Dynamics technique the bulk of the computa- 
tional cost is associated with the calculation of the far-field many-body interactions. However, 
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Figure 3: A cartoon of a dicolloidal particle. Radii of the two component spheres forming the 
dicolloid are a and b, while the distance between the centers of the two component spheres is 
c. In the present work only symmetric dicolloids with a = b are considered. 

particularly in concentrated suspensions, it is the lubrication interactions that dominate the 
overall particle-particle interactions. These considerations lead us to the development of the 
Fast Lubrication Dynamics technique (FLD), which is described below in the case of sus- 
pensions of spheres; its extension to anisotropic dicolloidal particles is presented later in Sec. 
(2.1.1). 

In the FLD approach, the overall resistance tensor is approximated as^*^'^^ 

R = Ro + R5, (6) 

where R5 is based on the near- field lubrication interactions summed pairwise^^'^^, while a 
novel isotropic term Rq (diagonal tensor, function of volume fraction alone) is selected to 
match the mean particle mobility (equivalently the short time self-diffusivity) from the more 
detailed Stokesian Dynamics technique for spheres Somewhat related approaches have 
been described before in the literature by Ball and Melrose and by Banchio and Brady 
The FLD approach has been shown to produce results comparable to the Stokesian Dynamics 
technique, while requiring only a fraction of the latter's cost; see Kumar and Kumar and 
Higdon^'^ for details; further validation is provided in the recent work of Schunk et al. We 
discuss next the extension of the FLD technique for modeling hydrodynamic interactions in 
dicolloidal particle suspensions. 

2.1.1 FLD technique for dicolloidal Particles 

In the FLD technique for dicolloidal particles (Fig. 3 shows a generic dicolloidal particle), 
the overall resistance tensor R is expressed by an equation similar to Eq. (6). For simplicity, 
in the present effort, we take the value of Rq in a suspension of dicolloids to be the same 
as in a suspension of spheres at the same volume fraction, irrespective of the exact shape 
of the particle under consideration. For particles with small degree of anisotropy this is a 
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Figure 4: A cartoon showing two pairwise interactions between spherical components forming 
the two dicolloids. 

fair assumption, though some error in the value of Rq can be expected at larger degree of 
anisotropy. Note that the largest aspect ratio particle investigated in the current work is 
AR = 1.5. Additionally, since our interest is only in concentrated suspensions with volume 
fraction (p > 0.3, we estimate the contribution from this error to be small, as the hydrodynamic 
interactions will be dominated by the strong lubrication forces at these high volume fractions. 

There are some additional considerations with the present choice of Rq, which deserve a 
brief comment. The tensor Rq in some ways represents the effective one-body resistance in 
the suspension. Our choice of Rq as an isotropic diagonal tensor implies that a dicolloidal 
particle, in the absence of interactions with other particles, will behave just like a sphere. 
For example, if we study an isolated dicolloidal particle with the current implementation of 
the FLD algorithm, then the dicolloid will not exhibit preferred flow alignment as might be 
expected for an anisotropic particle However, this is not a great concern here as the 
particle dynamics, including its orientation behavior, will be dominated by the lubrication 
interactions at the volume fractions considered here. 

Next, we discuss the computation of the near field lubrication term R^ in suspensions of 
dicolloidal particles. As shown in Fig. (3), a dicolloidal particle consists of two fused spheres 
(only symmetric dicolloids with a = b are considered here). Given this fact, the lubrication 
interaction between two dicolloids can be modeled as a pairwise sum of one or more interactions 
between the constituent spheres forming the two dicolloids; this is schematically shown in Fig. 
(4) and discussed in detail elsewhere Further note that up to four spherical node pairs can 
be formed between two dicolloidal particles, . In our algorithm, we do not necessarily consider 
all the four pairs. The selection of a spherical node pair for the lubrication calculation depends 
on two factors. First involves the satisfaction of a cut-off separation requirement - only those 
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spherical nodes pairs are considered in the lubrication calculation which have a center to center 
separation less than 2.5a, where a is the radius of the spherical node. Second factor involves the 
concept of inter-visibility - only those spherical node pairs are considered for the lubrication 
calculation which are inter- visible. Mathematically, if the line joining the centers of the two 
spherical nodes passes through each of their surfaces, then that pair has inter-visibility and 
is considered for the lubrication calculation; otherwise it is rejected. This scheme essentially 
prevents double counting which can happen otherwise. We remark that similar procedures 
are used extensively in the literature for computing interparticle interactions between patchy 
particles - a popular example is the Kern-Frenkel model 

2.2 Brownian Force (F^) 

Brownian forces and torques can be obtained from the fluctuation dissipation theorem and 
the equipartition of energy which dictates 

(F^) = (7a) 

(F^F^) = 2kT Kpu/^t (7b) 

where () denotes an ensemble average, R_f;7 denotes the component of the resistance ten- 
sor representing generalized force-velocity coupling, while At, k and T denote the time step 
employed in the algorithm, the Boltzmann constant, and the system temperature respec- 
tively^^'^^'^'^. Brownian forces/torques satisfying the above equations are easily obtained by 
employing uncorrelated random numbers with zero mean and unit variance. Briefly, we begin 
by writing the Brownian force/torque vector F^ as 



F^ = y^(Aa + B/3), (8) 

where a and /3 are 6Np vectors whose elements are uncorrelated random numbers, say ji, 
satisfying the aforementioned properties, namely 

= 0, (9a) 

{7i7j)=Sij, (9b) 

where 6ij is the Kronecker delta function. With the random numbers satisfying the above 
properties, we only need to satisfy the following two equations to obtain the agreement of Eq. 
(8) with Eq. (7) 



BB^ = R^^ 



(10) 

■s 
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The above equations require one to compute the square root of the matrices Rq'^ and Rf ^• 
In the present algorithm these are easily obtained as Rq^^ is a diagonal matrix, while contri- 
butions from JHf^ matrix can be summed pairwise. Details can be found in^^'^^. 

2.3 Repulsive Force 

During time stepping, due to integration errors, particle overlaps can occur. To overcome this, 
we employ a short range repulsive force between the spherical node pairs as follows^: 



Cp{^-^) {V-Iv'fd if 5<<5™„ 

(11) 

if 6 > 6min 



where 5 is the gap between the two spherical nodes, tj = 1 — 5/5min, d is the unit vector 
pointing from one spherical node to the other, and Cp is a parameter described below. In 
this work, the range of the repulsive force 6min was fixed at W~^a. For gaps smaller than a 
specified numerical tolerance, Snum, both the hydrodynamic force and interparticle force are 
capped by evaluating at 6num- This tolerance for all simulations presented here was set at 
^num = 10~^a. The interparticle force between a pair of dicolloids was computed as the sum 
of pairwise forces between the spherical nodes comprising the particles. The above form of 
the repulsive force gives a near hard sphere behavior with nonzero repulsive interaction only 
when the gap between a spherical node pair becomes less than 10~^a. 



2.4 Solution procedure 

The goal of the overall solution procedure is to compute the unknown particle velocities and 
angular velocities for a given set of Brownian forces, interparticle forces, and the imposed 
bulk (shear) flow. The unknown particle velocities/angular velocities can be obtained from 
the solution of the governing equation (3) written as: 

Rfc/ • (U°° - [/) = -(F^ + F^) -Rj7i5 •E°°. (12) 

The above equation is obtained from Eq. (3) by replacing F^ with its equivalent expressed 
in terms of the force-velocity (Ri?;/) and force-rate of strain (Rfe) couplings. In the present 
effort, we solve Eq. (12) iteratively with the conjugate gradient algorithm which is appli- 
cable here owing to the symmetry and positive definiteness of the resistance coupling Hpu- 
Once the particle velocities and angular velocities have been computed, the particle positions 
and orientations are evolved in time with a third order explicit Runge-Kutta method. The 
choices of the time-step At employed in the time integration step as well as the total time of 
integration t are discussed below. It is relevant to note here that all simulations in the present 
effort were initiated with a random configuration generated by a Monte Carlo method. A 
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Figure 5: Schematic of the system geometry: a cubic simulation box (side L) with Np = 1000 
particles subjected to a simple shear flow is considered. Periodic boundary conditions are 
employed in all three directions. 



Monte Carlo simulation for dicolloidal particles involved a random displacement to both the 
particle's orientation as well as its center of mass at each step, and this was continued until 
an apparent steady state was attained. 

2.5 Simulation Parameters 

In this section, we discuss some of important physical as well as numerical parameters in the 
current study. Perhaps the most important parameter in the present effort is the shear rate 7, 
which is expressed in terms of the non-dimensional Peclet number {Pe). Note that the Peclet 
number gives the ratio of the diffusive and the convective time scales, and is defined in the 
present effort as: 

(13) 

where a is the radius of the spherical node forming the dicolloid, and Dq = kT/Girfia is the 
Stokes-Einstein diffusivity based on the radius of the spherical node. Alternative definitions 
of Pe based on length scales other than a is possible, though that will not change the value 
of Pe appreciably given the small degree of anisotropy of dicolloids. In the current effort, Pe 
was varied between < Pe < 1000. Another important parameter in the current study is 
the particle volume fraction (p, which was varied in the range 0.3 < (j) < 0.5. Lastly, in the 
context of the present work, two particle shape parameters are required for defining a dicolloid. 
The two shape parameters are the ratio of the radii of the component spheres, b/a, and the 
ratio of the center to center distance and one of the radius, c/a; see Fig. (3) illustrating a 
generic dicolloidal particle. In the present study, we have investigated two different dicolloidal 
particles named as homonuclear and fused-dumbbell particles (Fig. 2). The shape parameters 
(b/a,c/a) for these two particle shapes are (1,0.2) for homonuclear particles and (1,1) for fused- 
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dumbbells. The geometric aspect ratio for these particles can be obtained as AR = (2a+c)/2a, 
which yields AR = 1.1 for homonuclear particles and AR = 1.5 for fused-dumbbell particles. 

There are several other parameters associated with the numerical solution procedure. First 
is the choice of the simulation box, which is taken to be cubic here with periodic boundary 
conditions employed in all three directions. A schematic of the system geometry is illustrated 
in Fig. (5). The number of particles in the simulation box was set to Np = 1000. The 
next parameter is the non-dimensional time-step At*, which varied in the range 10~^ — 10~^; 
typically, a lower At* is required at high volume fractions or low shear rates. We note that 
the non-dimensionalization of the time-step depends on the shear rate (Pe). At low shear 
rates {Pe < 1) , At is non-dimensionalized by the diffusive time scale, i.e.. At* = Do/a'^At, 
while at higher shear rates {Pe > 1), At is non-dimensionalized by the convective time scale, 
i.e.. At* = 'J At. The total non-dimensional integration time t* was t* = 500, with t being 
non-dimensionalized in the same way as At. This total integration time was found to be 
sufficient to obtain good statistical averages of various properties in the apparent steady 
state - in most cases, an apparent steady state was evident for times greater than t* > 
150 as evidenced by various of the suspension properties, like the orientation distribution 
function. The last parameter of interest is the repulsive force prefactor Cp, which in the non- 
dimensional form is set to C* = 10. At low shear rates {Pe < 1), Cp is non-dimensionalized 
as: C* = Cp/ {Sit ^a^ / Dq) , while at higher shear rates {Pe > 1), Cp is non-dimensionalization 
as: Cp = Cp/(87r/ia^7). The specification of various of the above parameters concludes the 
formulation section. We now turn our attention to the results from the numerical simulations, 
which are presented in the sections below. 

3 Results and Discussion 
3.1 Microstructure 

In this section, we present results for the positional (center-of-mass) microstructure in suspen- 
sions of homonuclear and fused-dumbbell particles. Results are presented for volume fractions 
between cf) = 0.3 and (p = 0.5 and for a wide range of Peclet numbers varying between Pe = 
and Pe = 1000. Figure (6) shows some illustrative simulation snapshots at Pe = 10 for both 
the homonuclear and the fused-dumbbell particles in (p = 0.3 and (f) = 0.5 suspensions. These 
snapshots are in the front view, i.e., in the velocity-gradient {x — y) plane. At the lower 
volume fraction of (/> = 0.3, there are no signs of ordering for both the particle shapes. At the 
higher volume fraction of (p = 0.5, several regions with a string like structure aligned in the 
flow direction are apparent in the homonuclear particle suspension, though the fused-dumbbell 
particle suspension appears mostly disordered even at this volume fraction. 

For a more quantitative characterization of the microstructure, we will employ the com- 
monly used measure of pair distribution function (PDF) denoted g{x,y,z). Note that the 
PDF gives the conditional probability of finding a particle at a given location {x, y, z) pro- 
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4) = 0.3 4) = 0.5 




(b) fused-dummbell 

Figure 6: Simulation snapshots in the front view {x — y plane) at a shear rate of Pe = 10 and 
for two different volume fractions {(j) = 0.3 and (j) = 0.5) for (a) homonuclear particles and (b) 
fused-dumbbell particles. 



vided a particle is already present at the origin. Additionally, the PDFs are normalized such 
that g{x, y,z) = 1 in the absence of any correlations between particle positions. Given the full 
three dimensional PDF g{x,y,z), it is usually instructive to consider its projections in various 
two dimensional planes - the projection in the front view [x — y plane) denoted g{x,y) will 
be of particular interest. Formally, g{x, y) is defined as follows: 

1 

5(2;,?/) = — y g{x,y,z)dz. (14) 

The PDF in the front view for homonuclear particles in (p = 0.3 and (f) = 0.5 suspensions 
are shown in Fig. (7) at several Peclet numbers between Pe = and Pe = 1000; figure 
(8) illustrates the same for fused-dumbbell particle suspensions. It is appropriate to remark 
here that the color scheme used in the PDF plot is chosen to emphasize small differences 
from the value of one, and, as such, different patterns in these plots do not necessarily imply 
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Figure 8: Pair distribution function in the front view {x — y plane) for fused-dumbbell particles 
at several Pe and two different volume fractions {(j) = 0.3 and (p = 0.5). 
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magnitudes greatly different from one. Returning to the above two figures (Figs. 7 and 8), 
we note that at equihbrium (Pe = 0), a disordered microstructure is apparent in suspensions 
of both the homonuclear and fused-dumbbell particles at (p = 0.3; in these plots, only the 
first coordination shell distinctive of an amorphous fluid like structure is apparent. At the 
higher volume fraction of = 0.5, the structure stays disordered at Pe = 0, though in these 
cases the second (for fused-dumbbell) or third (for homonuclear) coordination shells are also 
apparent. A disordered structure at equilibrium at these volume fractions is consistent with 
the Monte Carlo simulations of Vega et al. . As per the previous reference, a suspension of 
homonuclear particles are in a fluid like state below a volume fraction of = 0.4945 (called 
freezing point) and a crystalline solid state above (j) = 0.543 (called melting point); between 
these two volume fractions the two phases coexist. In fused-dumbbell particle suspensions, 
the freezing and melting transitions were found to occur at approximately 4> = 0.6013 and 
(j) = 0.6463. Typically, in simulations, the ordered phase is observed at volume fractions 
above the melting transition Therefore, a fluid like state is expected at equilibrium in 
suspensions of homonuclear as well as fused-dumbbell particles for volume fractions up to 
(p = 0.5 considered here. It is interesting to remark here that the ordered phase at the freezing 
transition is orientationally disordered (plastic crystal) face centered cubic in suspensions of 
homonuclear particles, while it is orientationally ordered base centered monoclinic in fused- 
dumbbell particle suspensions'^^. It is also pertinent to remark here that a nematic phase is 
not expected in dicolloidal particle suspensions. This is due to the fact that the degree of 
anisotropy of dicolloids is small, while the nematic phase is stable only at higher aspect ratios. 

As an example, for spherocylinders, the nematic phase becomes stable at an aspect ratio of 
4.733^ 

We next discuss the microstructure in the (p = 0-3 suspensions of homonuclear and fused- 
dumbell particles at finite rates of shear. At very low shear rates (Pe = 0.01), the microstruc- 
ture is found to be very similar to the equilibrium microstructure in both the homonuclear 
and the fused-dumbbell particle suspensions (not shown). As the shear rate is further in- 
creased {Pe = 1 in the first column of Figs. 7 and 8), the microstructure becomes distorted 
- there is an enhancement in pair probability in the compressional quadrants, accompanied 
by its depletion in the extensional quadrants. This fore- aft asymmetry in the PDF results 
due to the presence of Brownian forces as discussed in Brady and Morris and Foss and 
Brady Lastly, we note that the microstructure remains disordered at the highest shear 
rates investigated here {Pe = 1000) in the (p = 0.3 suspensions (Figs. 7 and 8). 

At the higher volume fraction of = 0.5, some interesting features emerge in the PDFs of 
both the particle shapes (Figs. 7 and 8). At low shear rates {Pe = 0.01), the microstructure is 
similar to the equilibrium microstructure for both the particle shapes (not shown). However, 
at Pe = 10, evidence of weak string like ordering is apparent in suspensions of homonuclear 
particles (horizontal bands in Fig. 7), which is sustained at higher shear rates, though the 
ordering become progressively weaker (e.g., compare PDF at Pe = 10 with PDF at Pe = 
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1000). In fused-dumbbell particle suspensions, the microstructure shows similar hints of string 
like ordering at Pe = 10, though the strength of ordering is considerably weaker than in 
suspensions of homonuclear particles (Fig. 8). At further higher shear rates in fused-dumbbell 
particle suspensions (e.g., Pe = 1000), the microstructure appears mostly disordered. 

The results presented in the current study are consistent with previous experimental and 
computational efforts. In experimental studies on near hard sphere suspensions, Ackerson^^ 
noted that suspensions which are disordered (liquid-like) at rest stays disordered under shear. 
This is consistent with the disordered structure observed in = 0.3 suspensions of homonuclear 
particles under shear, and also in fused-dumbbell particle suspensions for volume fractions 
up to (/> = 0.5 considered here. Ackerson^^ also noted that at volume fractions near the 
freezing transition, string like structures begin to coexist with amorphous fluid like regions 
at intermediate rates of shear - this is similar to the weak string like ordering observed in 
(p = 0.5 homonuclear particle suspension. Similarly, in Stokesian Dynamics simulations on 
non-Brownian dicolloidal particle suspensions, which is equivalent to the high Pe limit, Kumar 
and Higdon^ also observed weak string like ordering in suspensions of homonuclear particles 
at (p = 0.5. In independent Stokesian dynamics simulations of Sierou and Brady again for 
non-Brownian particles, very strong string like ordering was observed in suspensions of spheres 
at (j) = 0.52, which is just above the freezing point of = 0.494 in hard sphere suspensions. It 
is instructive to note that in suspensions that are ordered at rest, like in sufficiently charged 
particle suspensions, very distinct and rich set of microstructural transitions are observed as 
a function of shear rate; see, e.g., Chen et al.^^. Chow and Zukoski^'^. This regime, however, 
is not relevant to our study here. Finally, there have been very limited experimental studies 
on the microstructure of dicolloidal particle suspensions under steady shear. The only study 
we are aware of is that by Mock and Zukoski'^^, who studied charge- stabilized dicolloidal 
particle suspensions under shear at a volume fraction of (/> = 0.42. The homonuclear particle 
equivalent investigated in this work was found to possess a polycrystalline structure at rest, a 
sliding layer structure at intermediate rates of shear, and an amorphous structure at further 
higher shear rates. A direct comparison of this work with the present eff'ort, though, is not 
possible due to the presence of strong electrostatic repulsion in the former study. 

3.2 Orientation behavior 

We next present the results on the orientation behavior in homonuclear and fused-dumbbell 
particle suspensions as a function of Pe for volume fractions between (f> = 0.3 and 4> = 0.5. 
The statistical distribution of the particle orientation will be characterized by the orientation 
distribution function (ODF) denoted P{9,ip), where 9 and (p are spherical coordinates as 
illustrated in Fig. (2); also see Fig. (10) for illustrations of the angle 99. In words, 9 is the 
polar angle between the director n of the particle and the z-axis, while (p is the azimuthal 
angle between its projection in the x — y plane and the y-axis measured clockwise. Based on 
the ODF P{9, ip), the differential probability of finding a particle with an orientation between 



18 



(9, 9+d9) and {(f,(f+d(p), dN, can be expressed as follows: 



dN = ^ P{e, (p) sin e de dip. (15) 
47r 



Additionally, we note that the P{9,(p) satisfies the following normalization constraint: 

f7r/2 /■7r/2 
t/2 



/ P{e,ip) sine de dip = 1, (16) 
-it/2 Jo 



where we have restricted the limits of integration in e to [0, 7r/2] and in ip to [— 7r/2, 7r/2], 
which is half of the corresponding overall ranges; this is sufficient here owing to the symmetry 
of the ODF P{e,ip). We next introduce one dimensional ODFs, P{e) and P{ip), which are 
defined as follows: 

p{e)= r Pie, ip) dip, (17a) 

J-Tr/2 
pw/2 

P{ip) = iT I P{e,ip) sinede. (17b) 
Jo 

We will use P{e) and P{(p) extensively in lieu of the full ODF P(9,(p). It is important to 
emphasize here that the definitions of P{e) and P{(p) are such that a uniform distribution in 
the orientation space will result in P(e) = 1 and P{ip) = 1. For later use, we also define a 
suspension averaged e and ip, denoted and ipm, as follows: 

1-77/2 

em= 9P{e) sinede, (ISa) 
Jo 

1 

iPm = - \ip\P{ip)dip. (18b) 

7!" J-n/2 

Note the absolute value term '|(^|' in the definition of ipm above. A simple interpretation of 
em is that it gives the mean inclination of the particle director with the z axis. Similarly, ip^, 
gives the mean inclination of the projection of the particle director in the x — y plane with 
the y axis. As should be obvious from the above definitions, e^ and ipm have the following 
bounds: < ^ and < ipm < tt/2. In a suspension with a completely random 
orientation distribution, these two averages will be denoted by e^ and ip^- We then define 
Aem and Aip^ for an arbitrary suspension as follows: 

Aem = em- C (19a) 
Aipm = 'Pm-^m- (19b) 

When Aem is negative, it will imply an enhanced vorticity alignment in comparison to an 
isotropic distribution, while positive Aem will imply an enhanced alignment in the velocity- 
gradient plane (i.e., away from the vorticity axis) in comparison to an isotropic distribution. A 
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positive A(^m will imply an enhanced alignment with the velocity axis in the velocity-gradient 
plane, while a negative A(/7m will imply an enhanced alignment with the gradient axis in the 
velocity-gradient plane, both, of course, are relative to an isotropic distribution. 

In addition to the above measures, it will be useful to employ direct measures of the 
particle orientation relative to the three coordinate axes x, y, and z. In the present work, 
this measure is provided by the orientation order parameter Su defined relative to the axis i, 
where i could be x, y, or z. More specifically, the orientational order parameter Su is defined 
as the average over all particles of 

3 1 

Sii = -niUi - -. (20) 

Some representative values of Su are as follows. For a random orientation distribution Su = 0; 
for perfect alignment of particles with axis Su = 1.0; and for particles aligned in the plane 
perpendicular to axis i, Su = —0.5. It is evident from the above definition that an increased 
degree of particle alignment with an axis i will result in a higher value of Su- 

Having introduced the various measures for characterizing the orientation behavior, we 
turn our attention to the results for these measures. We will first present the orientation 
distribution in the 6 coordinate as a function of Pe, volume fraction and particle shapes 
in Sec. (3.2.1). Following this, we will present the orientation distribution in ip coordinate 
in Sec. (3.2.2). Lastly, we will present the results for the orientation order parameters in 
Sec. (3.2.3). 

3.2.1 Orientation distribution in 6 coordinate 

Figures (9a) and (9b) show the orientation distribution function P{9) at several Peclet numbers 
in homonuclear particle suspensions with = 0.3 and (j) = 0.5, respectively. The same two 
plots for suspensions of fused-dumbbell particles are shown in Figs. (9d) and (9e), respectively. 
The change in the suspension averaged 9 over an isotropic distribution, /S.9m, is shown in 
Fig. (9c) for homonuclear particle suspensions, while it is shown in Fig. (9f) for fused-dumbbell 
particle suspensions. 

We first consider the orientation behavior in suspensions of homonuclear particles. At low 
shear rates {Pe = 0.01) in the ^ = 0.3 suspension, P{9) has a value close to one, implying a 
random distribution in 9 at these shear rates (Fig. 9a). As the shear rate is further increased, 
an enhancement in probability at low 9 values is observed, indicating an increased vorticity 
alignment; this enhancement is particularly strong at Pe = 1000 (Fig. 9a). We see a very 
similar trend at the higher volume fraction of = 0.5 (Fig. 9b), namely, a random distribution 
in 9 at low shear rates {Pe = 0.01), which shifts towards lower 9 values, i.e., increased vorticity 
alignment, at higher shear rates (e.g. Pe = 1000). These trends in the orientation behavior 
with shear rate are also evident in the plot of in Fig. (9c). At low shear rates, A^^ is 
nearly zero indicating a random orientation distribution, while it becomes progressively more 
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Figure 9: P{0) and A0m for (a), (b), and (c) homonuclear particles and (d), (e) and (f) 
fused-dumbbell particles. 



negative with increasing shear rates, implying an increased vorticity alignment at higher shear 
rates. It is also evident in the plot that A6m is more negative at higher volume fractions, 
implying a greater vorticity alignment at higher volume fractions in homonuclear particle 
suspensions. 

We next discuss the orientation results in suspensions of fused-dumbbell particles. Focusing 
on the plot of P{0) in cp = 0.3 fused-dumbbell particle suspension (Fig. 9d), it is obvious that 
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the distribution at low shear rates {Pe = 0.01) is close to random. With increasing shear rates, 
however, P{0) changes from a random distribution to a distribution skewed towards higher 9 
values, indicating an enhanced alignment in the velocity-gradient plane (Fig. 9d). This trend 
is clearly opposite to the trend observed in homonuclear particle suspensions. At the higher 
volume fraction [<j) = 0.5), the behavior is more complex (Fig. 9e). In this case, at intermediate 
shear rates {Pe = 1), an enhanced alignment in the velocity-gradient plane (i.e., at 9 values) 
is observed. With further increase in shear rate, however, the enhanced particle alignment in 
the velocity-gradient plane is attenuated (e.g., at Pe = 1000), and the distribution gradually 
shifts towards vorticity alignment, though it is not as complete or as strong as in homonuclear 
particle suspensions. We confirm these observations by examining the plot of in Fig. (9f). 
It is evident from the figure that, at low volume fractions {(j) = 0.3), IS.9m increases nearly 
monotonically with increasing shear rate, indicating an increasing alignment in the velocity- 
gradient plane. However, at (/> = 0.5, initially increases with increasing shear rate till 
Pe = 1, and thereafter it decreases and becomes negative with increasing shear rates before 
settling at a near zero value. 

In a recent study ^, we investigated the orientation distribution in the 9 coordinate in 
sheared non-Brownian suspensions of dicolloids. More precisely, the orientation distribution 
into Jeffery's orbits was investigated, which approximately corresponds to a distribution in 
the 9 coordinate for small aspect ratio particles like dicolloids. This study indicates that in 
low volume fraction suspensions, both the homonuclear and fused-dumbbell particles exhibit 
increased particle alignment in the velocity-gradient plane. The increase in the velocity- 
gradient alignment in suspensions of fused-dumbbell particles was found to be considerably 
higher than in the suspensions of homonuclear particles. At higher volume fractions, both the 
particle shapes were found to exhibit a shift towards enhanced vorticity alignment, though 
the degree of vorticity alignment was much weaker in suspensions of fused-dumbbells than in 
suspensions of homonuclear particles. The origins of the contrasting behavior at low and high 
volume fractions was attributed to two different mechanisms at play in these two parameter 
regimes. At low volume fractions, on an average, it was observed that a particle aligned with 
the vorticity axis underwent comparatively stronger fluctuation in its orientation during the 
course of a collision with another particle. This behavior leads to non-uniform fluctuations 
in the orientation space with a stronger fluctuation observed for vorticity aligned particles. 
Additionally, the non-uniformity of fluctuations was found to be stronger and distributed over 
a wider region in the orientation space for the more anisotropic fused-dumbbell particles than 
the homonuclear particles. As discussed in Shaqfeh and Koch^'', these non-uniform fluctua- 
tions will result in a drift away from regions of higher fluctuations towards regions of lower 
fluctuations. Predictions based on this mechanism easily resolves two of the key observations 
at low volume fractions, namely, an increased alignment in the velocity-gradient plane in both 
the homonuclear and the fused-dumbbell particle suspensions, and a comparatively greater 
alignment in the velocity-gradient plane in fused-dumbbell particle suspensions. At higher 
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volume fractions, the fluctuations were found to be nearly uniform in the orientation space 
and hence cannot lead to a non-uniform orientation distribution. The high volume fraction 
behavior was instead found to be driven by the negative second normal stress difference ex- 
hibited by these suspensions at high volume fractions (note that at low volume fractions, the 
normal stress differences are very weak). A detailed discussion of this mechanism is beyond 
the scope of the current work and the interested reader is referred to the original publication^. 
The greater vorticity alignment of the homonuclear particles than that of the fused-dumbbell 
particles at high volume fractions was shown to result from the lower degree of anisotropy of 
homonuclear particles. In essence this implies that the homonuclear particles have a greater 
rotational "mobility" than fused-dumbbell particles in concentrated suspensions, and hence 
the second normal stress difference is a more effective driving force in the former case. 

With the aid of the above discussion, we can provide some mechanistic insights into the 
observed orientation behavior in the current study. Firstly, in all cases, a random orientation at 
low shear rates is expected due to the dominant effect of the Brownian forces in this regime. 
The high shear rate regime is well explained by the non-Brownian limit discussed above. 
The non-Brownian limit consistently explains the enhanced vorticity alignment observed in 
suspensions of homonuclear particles at high shear rates and at all volume fractions considered 
here. It also consistently explains the enhanced velocity-gradient plane alignment in low 
volume fraction suspensions of fused-dumbbell particles, and the (weak) vorticity alignment 
in higher volume fraction suspensions of fused-dumbbell particles, both in the high shear 
rate limit. The only regime that remains to be resolved is the intermediate shear rate regime. 
Intuitively, one may expect that the behavior at intermediate shear rates will be "intermediate" 
of the high and the low shear rate behaviors. This is certainly found to be true, but there is 
one exception - the intermediate shear rate behavior in the (p = 0.5 fused-dumbbell particle 
suspension (Fig. 9f). In this case a random distribution is observed at low shear rates, while 
a weak vorticity alignment is observed at higher shear rate. The intermediate shear rate, 
however, exhibits enhanced velocity-gradient alignment. The origins of this behavior is not 
entirely clear, though it is possible that the "non-uniform fluctuation" mechanism plays a role 
at intermediate shear rates in high volume fraction suspensions as well (at very low shear 
rates Brownian motion will dominate any other effect). Further studies are clearly necessary 
to confirm this mechanism. 

3.2.2 Orientation distribution in (p coordinate 

The orientation distribution in the <p coordinate will be characterized by the ODF P{<p) as 
well as by Aipm, the change in the suspension averaged ip in comparison to an isotropic 
distribution. For homonuclear particles, Fig. (11a) shows P{ip) in (p = 0.3 suspension and 
Fig. (lib) shows the same in = 0.5 suspension; the plot for Aipm is shown in Fig. (11c). For 
fused-dumbbell particles. Fig. (lid) shows P{ip) in = 0.3 suspension and Fig. (lie) shows 
the same in = 0.5 suspension; the plot for Aiprn is shown in Fig. (llf). 
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Figure 10: Examples of particles aligned in the compressional and the extensional quadrants. 

We begin our discussion by examining the ODF -P((/5) in homonuclear particle suspensions 
at (j) = 0.3 (Fig. 11a). At low shear rates {Pe = 0.01), a random distribution is obvious 
as P{(p) ~ 1. At intermediate rates of shear (Pe = 1) for the same suspension, we find a 
decrease in probability corresponding to orientations in the compressional quadrant (—0.5 < 
if/n < 0) accompanied by an increase in probability corresponding to orientations in the 
extensional quadrant (—0.5 < ip/n < 0). Examples of particles aligned in the compressional 
and extensional quadrants are illustrated in Fig. (10). At higher shear rates (e.g., Pe = 1000 
in Fig. 11a), the height of probability maxima as well as the depth of probability minima 
increases. Additionally, the probability minima is found to shift closer to the gradient axis 
{ip = 0) and the probability maxima closer to the velocity axis {ip = tt/2) with increasing 
shear rates (e.g., compare Pe = 1 and Pe = 1000 in Fig. 11a). Turning our attention to 
suspensions of fused-dumbbell particles at the same volume fraction (</) = 0.3), we observe the 
same trends in the ODF P{(p) as in the suspension of homonuclear particles (compare Fig. 
11a and lid). 

At the higher volume fraction of = 0.5, the ODF P{(p) shows additional features in 
both the homonuclear and the fused-dumbbell particle suspensions (Figs, lib and lie). In 
the homonuclear particle suspension, an expected random distribution is evident at Pe = 0.01 
(Fig. lib). At intermediate shear rates {Pe = 1), similar to the (p = 0.3 suspension, an 
enhancement in probability corresponding to orientations in the extensional quadrant and 
a reduction in probability corresponding to orientations in the compressional quadrant is 
observed (Fig. lib). However, with further increase in shear rate {Pe = 10), the height of the 
probability maxima as well as the depth of the probability minima decreases (Fig. lib). This 
is in contrast to the behavior at low volume fractions discussed above. This trend at Pe = 10 
is also evident in fused-dumbbell particle suspensions at this volume fraction (Fig. lie). At 
further higher shear rates (e.g., Pe = 1000), the height of the maxima as well as the depth 
of the minima generally increases in suspensions of the homonuclear particles (Fig. lib), but 
decreases slightly in suspensions of fused-dumbbell particles (Fig. lie). 
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Figure 11: P{(p) and A(^m for (a), (b), and (c) homonuclear particles and (d), (e) and (f) 
fused-dumbbell particles. 



The trends in the ODF P{(p) discussed above are confirmed in the Aipm plots shown in 
Fig. (11c) for homonuclear particles and in Fig. (llf) for fused-dumbbell particles. For both 
homonuclear and fused-dumbbell particles, there is a general tendency of increasing Aipm with 
increasing volume fractions. Further, for both particle shapes, Aipm increases with increasing 
shear rates at low volume fractions (e.g., (p = 0.3). At high volume fractions, however, Aipm 
does not necessarily increase with increasing shear rates. This is particularly evident in the 
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Figure 12: A cartoon showing a slender body subjected to an extensional flow. Also shown is 
a small sphere "colliding" with the slender body. 

(p = 0.5 homonuclear particle suspensions, where A(/9m is found to decrease in a small shear rate 
window around Pe = 10 (Fig. 11c). In fused-dumbbell suspensions at high volume fractions, 
Ac^m is found to decrease at all shear rates above Pe > 1 (Fig. llf). Another feature apparent 
in the plots of Aiprn for the two particle shapes is that its magnitude is considerably larger 
in suspensions of fused-dumbell particles than in suspensions of homonuclear particles. It is 
worth recalling at this point that Aipm has important implications for results on the flow and 
the gradient alignment of the particle, e.g., a higher positive Aipm will imply a greater degree 
of flow alignment at the expense of gradient alignment in the velocity-gradient plane. These 
aspects are discussed further in Sec. (3.2.3). 

Some of the trends in the ODF P{(p) as a function of shear rate have striking similarities 
with the corresponding trends in suspensions of Brownian axisymmetric particles at infinite 
dilution. We discussed the orientation behavior in such systems in the introduction and we 
summarize the key observations here. At equilibrium, the orientation distribution is expectedly 
isotropic. As the suspension is sheared at very low shear rates, to leading order, the vorticity 
component of the shear flow has no effect on the orientation distribution due to the isotropic 
distribution at equilibrium. The effect of the extensional component of the flow is to rotate the 
particle away from the compressional axis towards an alignment with the extensional axis. As 
a consequence, the probability near the compressional axis is reduced, while the probability 
near the extensional axis is enhanced The higher order approximations result from the 
vorticity as well as the extensional component of the flow^^. In particular, the effect of the 
vorticity component of the flow is to rotate the distorted orientation distribution clockwise, 
such the probability maxima in the ODF shifts closer to the velocity axis, while the probability 
minima shifts closer to the gradient axis. 

It is evident that the ODF P{ip) in the present effort exhibits some trends similar to that 
in suspensions at infinite dilution. However, the mechanism of the distortion of the ODF 
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upon the onset of shear is very different in the two cases. First of all, we remind the reader 
that the effective single particle resistance tensor Rq of dicolloids in FLD is diagonal, which is 
qualitatively similar to the resistance tensor of an isolated sphere. Hence a single dicolloid by 
itself cannot exhibit a preferred alignment with any axis upon being subjected to shear. So, 
the distortion of the ODF in the present case clearly results from particle-particle interactions, 
which are expected to dominate anyway in concentrated suspensions. The precise mechanism 
of this distortion is easily understood with the help of a highly idealized example presented 
below. 

Consider a single slender axisymmetric particle subjected to extensional flow, as shown 
in Fig. (12). This single particle by itself cannot rotate upon being subjected to extensional 
flow due to reasons discussed above. Now consider a second particle, which for simplicity is 
taken as a relatively small sphere in Fig. (12). The sphere, to first approximation, will move 
along the streamlines and consequently "collide" with the slender body causing it to rotate 
clockwise towards the extensional axis. It is easy to see that this will be case for most, if 
not all, starting positions of the second particle. Thus, in this idealized case, we see that 
the effect of particle collisions on the orientation dynamics of the slender body is much the 
same as that of slender body by itself in the extensional flow. By analogy of the orientation 
behavior at infinite dilution in Brenner , we may expect an enhanced orientation probability 
along the extensional axis and a reduced orientation probability along the compressional axis 
at low shear rates in the present case as well. This is obviously the case here (see Fig. 11). 
Furthermore, at higher shear rates (Pe), the distorted ODF is expected to rotate clockwise 
so as to shift the probability maxima closer to the velocity axis and the probability minima 
closer to the gradient axis. This again is very similar to the observations in the present effort. 
However, there are some additional features in the ODF P{(p) not resolved yet. At Pe = 10 
in the higher volume fraction suspensions of both the homonuclear and the fused-dumbbell 
particles, the depth of the minima as well as the height of the maxima of P{(p) decreases from 
the corresponding values at Pe = 1; see Figs. 11(b) and 11(e). It appears that this behavior 
results due to the formation of weak string like ordering around Pe = 10, particularly in 
suspensions of homonuclear particles. An obvious consequence of the presence of string like 
ordering will be to reduce the collision frequency between the particles, which, as discussed 
above, is a key driver of the distortion in the ODF. Yet another feature of the ODF P{ip) in 
fused-dumbbell particle suspensions at high volume fractions and high shear rates is that it 
continues to display a decrease in the height of probability maxima and the depth of probability 
minima with increasing Pe (Fig. lie). This cannot result from a string like ordering transition, 
as the microstructure is completely disordered in these suspensions at high shear rates (see 
Sec. 3.1). A likely origin of this behavior could be the hydrocluster mechanism proposed by 
Egres et al. - at high shear rates, large clusters of particles are typically formed, and these 
clusters of particles are expected to have a lower degree of anisotropy than the individual 
particles themselves. It is likely that the formation of more spherical clusters will result in a 
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Figure 13: Orientation order parameters as a function of shear rate for (a) and (b) homonuclear 
particles and (c) and (d) fused-dumbbell particles at two different volume fractions (</> = 0.3 
and (j) = 0.5). 



reduced distortion of the ODF. Lastly, we note that the trends in Aip^ follows directly from 
the corresponding trends in the ODF P{(p), and hence will not be explicitly discussed here. 

3.2.3 Orientation order parameters 

We discussed above the orientation distribution functions P{9) and P{f). These measures 
are particularly useful for gaining a mechanistic understanding of the orientation behavior. 
However, a more natural way of presenting the orientation results is to quantify the particle 
alignment with the three coordinate axes: velocity (x), gradient (y), and vorticity (z). With 
this goal, we present the results for the orientation order parameters Sxx, Syy, and Szz below. 

Figures 13(a) and 13(b) respectively show the orientation order parameters va. (p = 0.3 and 
(j) = 0.5 suspensions of homonuclear particles as a function of Pe. \n. (j) = 0.3 suspensions, at 
Pe = 0.01, all three orientation order parameters are nearly zero, consistent with a random 
orientation at this low shear rate. With increasing shear rates, Szz is found to increase, 
indicating an increased vorticity alignment with increasing shear rate as discussed in Sec. 
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3.2.1. In contrast, Sxx as well as Syy decreases with increasing shear rate; the decrease in Sxx 
is smaher than the decrease in Syy. The decrease in Sxx and Syy primarily result due to the 
increase in Szz^ which leaves a smaller fraction to be distributed between Sxx and Syy. The 
fact that Sxx decreases slower than Syy is a direct consequence of the orientation distribution 
in the </? coordinate - with increasing shear rate, flow alignment is favored at the expense 
of gradient alignment as discussed in Sec. 3.2.2. At the higher volume fraction ol (p = 0.5, 
the trends in orientation order parameters are very similar to the trends observed in i;^ = 0.3 
suspension - the order parameter S^z increases with shear rate, while Syy decreases with 
shear rate. The only difference is that, at = 0.5, Sxx does not decrease monotonically with 
increasing shear rate as is the case in the (f) = 0.3 suspension. As noted above, the instances in 
which Sxx increases with shear rate results due to the enhanced flow alignment of the particle 
in the velocity-gradient plane at the expense of the gradient alignment; see Sec. 3.2.2. 

We next consider the orientation order parameters for fused-dumbbell particles in = 0.3 
and (p = 0.5 suspensions, which are shown in Figs. 13(c) and 13(d), respectively. In the (j) = 0.3 
suspension, the order parameter Szz decreases with increasing shear rate, which is the direct 
consequence of increasing with increasing shear rate; see Sec. 3.2.1. The order parameter 
Sxx monotonically increases with increasing Pe. This results both due to the monotonically 
increasing with increasing Pe, implying an increasing velocity-gradient alignment with 
increasing Pe, as well as a (nearly) monotonically increasing A(/?m with increasing shear rate, 
indicating an increasing flow alignment in the velocity-gradient plane. The order parameter 
Syy initially decreases with increasing shear rate, which primarily results due to the increase 
in Sxx- However, with further increase in shear rate, Syy is found to increase, which can be 
attributed to the orientation drift away from the vorticity axis at higher shear rates, essen- 
tially leaving a greater fraction to be distributed in the flow and the gradient direction. At 
the higher volume fraction of = 0.5, the order parameter Szz^ consistent with the trends in 
A^mj initially decreases with increasing shear rates, but is then found to increase for Pe > 1 
and subsequently levels off at near zero values at very high shear rates. The order parameter 
Sxx shows a non-monotonic behavior with increasing shear rate: it initially increases with 
increasing shear rate reaching a maxima at Pe = 1; thereafter, it decreases with further in- 
crease in shear rate and eventually levels off at the highest shear rates considered here. In 
conjuction with the non-monotonic Sxx, the variation in the order parameter Syy with increas- 
ing shear rate is also non-monotonic, though in the opposite direction: Syy initially decreases 
with increasing shear rate, and thereafter it increases before it levels off at the high shear 
rate plateau. The non-monotonic variation in Sxx (or Syy) with shear rate primarily results 
from the corresponding non-monotonic variation of Aip^ with shear rate. As discussed in Sec. 
(3.2.2), the non-monotonic variation in Aipm perhaps results from two distinct mechanisms 
- first is the appearance of a weak string like ordering around Pe = 10, and second is the 
formation of hydroclusters at high shear rates . It is interesting to note that both of these 
mechanisms have the same apparent effect of causing a reduction in the flow alignment of the 
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particles. 



4 Conclusions 

In this study we investigated the orientation and microstructure in sheared Brownian suspen- 
sions of hard dicolloidal particles. Studies were conducted for a wide range of volume fractions 
(0.3 < 4> < 0.5) and shear rates (0 < Pe < 1000) in suspensions of homonuclear particles (as- 
pect ratio 1.1) and fused-dumbbell particles (aspect ratio 1.5). The microstructure was found 
to lack a long range order at all volume fractions, though evidence of a weak string like or- 
dering was apparent in suspensions of homonuclear particles at high volume fractions and at 
intermediate to high shear rates. The orientation distribution was found to be random at low 
shear rates. At the higher shear rates, homonuclear particles exhibited enhanced vorticity 
alignment at all volume fractions; fused-dumbbell particle suspensions, in contrast, exhibited 
reduced vorticity alignment at low volume fractions, and a negligible shift relative to the vor- 
ticity axis at higher volume fractions. The orientation distribution in the velocity-gradient 
plane indicated an increased flow alignment with increasing vohime fractions for both particle 
shapes; flow alignment also increased with increasing shear rates at low volume fractions, but 
not necessarily at higher volume fractions. These observations on the orientation distribution 
were consistently resolved in the light of previous studies on Brownian and non-Brownian 
suspensions. 
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